library(tidyverse); theme_set(theme_bw())
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(fs);
library(cowplot); theme_set(theme_cowplot(font_size=6));
##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:lubridate':
##
## stamp
library(trend)
# File paths
# roi_file <- "/labs/kingsley/ambenj/myosin_dups/analysis/transplant_pops/gasAcuv5_C4masked/06_samtools_coverage_roi/roi_coverage_mapq3.txt"
# wg_path <-"/labs/kingsley/ambenj/myosin_dups/analysis/transplant_pops/gasAcuv5_C4masked/06_samtools_coverage_wg"
# af_file <- "/labs/kingsley/ambenj/myosin_dups/analysis/transplant_pops/Veeramah_PoolSeq_AFs_corrected.filtered_chrXIX.txt"
# out_file <- "/labs/kingsley/ambenj/myosin_dups/analysis/transplant_pops/gasAcuv5_C4masked/poolseq_depth_figure.pdf"
# File paths
roi_file <- "/Volumes/lab_kingsley/ambenj/myosin_dups/analysis/transplant_pops/gasAcuv5_C4masked/06_samtools_coverage_roi/roi_coverage_mapq3.txt"
wg_path <-"/Volumes/lab_kingsley/ambenj/myosin_dups/analysis/transplant_pops/gasAcuv5_C4masked/06_samtools_coverage_wg"
af_file <- "/Volumes/lab_kingsley/ambenj/myosin_dups/analysis/transplant_pops/Veeramah_PoolSeq_AFs_corrected.filtered_chrXIX.txt"
out_file <- "/Volumes/lab_kingsley/ambenj/myosin_dups/analysis/transplant_pops/gasAcuv5_C4masked/poolseq_depth_figure.pdf"
out_file_af <-"/Volumes/lab_kingsley/ambenj/myosin_dups/analysis/transplant_pops/gasAcuv5_C4masked/poolseq_AF_figure.pdf"
# Read roi file
roi <- read_tsv(roi_file) %>%
mutate(bam = str_remove_all(bam, "05_mkdup_index/|.PE_ME_merged.sort.Mkdup.realign.BQrecal.realignGA5_C4masked.sort.merged.mkdup.bam|.PE_ME_merged.sort.realign.BQrecal.realignGA5_C4masked.sort.merged.mkdup.bam|.PE_ME_merged.sort.Mkdup.rename.realign.BQrecal.realignGA5_C4masked.sort.merged.mkdup.bam")) %>%
rename(samp = bam)
## Rows: 253 Columns: 12
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): desc, region, chr, bam
## dbl (8): startpos, endpos, numreads, covbases, coverage, meandepth, meanbase...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# separate(bam, into = c("pop", "year"), sep="-", convert = TRUE)
roi
## # A tibble: 253 × 12
## desc region chr startpos endpos numreads covbases coverage meandepth
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 NLRC5 chrXI… chrX… 2601498 2.63e6 15164 30076 95.3 67.9
## 2 MYH3C1 chrXI… chrX… 2637324 2.65e6 4770 13318 89.2 47.1
## 3 MYH3C2 chrXI… chrX… 2655758 2.67e6 3912 10519 96.8 54.3
## 4 MYH3C3 chrXI… chrX… 2667295 2.68e6 3602 10178 94.9 50.6
## 5 duplicatio… chrXI… chrX… 2665791 2.68e6 5039 16600 94.2 42.0
## 6 SYT19 chrXI… chrX… 2697955 2.73e6 11447 24674 75.5 50.9
## 7 CALB2A chrXI… chrX… 2745064 2.76e6 13620 14596 100 138.
## 8 HTRA1A chrVI… chrVI 14382418 1.44e7 25067 16913 98.7 217.
## 9 MYH3C3y chrY:… chrY 8892731 8.90e6 3282 11244 100 43.4
## 10 MYH3C2y chrY:… chrY 8904789 8.92e6 4912 10952 100 68.2
## # ℹ 243 more rows
## # ℹ 3 more variables: meanbaseq <dbl>, meanmapq <dbl>, samp <chr>
# Read whole genome files
files <- dir_ls(path = wg_path, glob = "*wg_coverage.txt")
# function to add file name to dataframe
read_and_record_filename <- function(filename){read_tsv(filename) %>%
mutate(filename = path_file(filename))
}
# gather files into one dataframe
wg <- map_df(files, read_and_record_filename) %>%
separate(filename, into=c("samp"), sep="\\.", extra="drop")
wg
## # A tibble: 552 × 10
## chr startpos endpos numreads covbases coverage meandepth meanbaseq meanmapq
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 chrI 1 2.96e7 37727036 29260296 98.8 190. 22.4 58.1
## 2 chrII 1 2.37e7 27651994 23455798 99.0 174. 22.5 58.6
## 3 chrI… 1 1.78e7 23018107 17487754 98.5 194. 22.5 57.8
## 4 chrIV 1 3.42e7 41943758 33526396 98.1 183. 22.5 57.9
## 5 chrIX 1 2.08e7 26353055 20426020 98.0 189. 22.5 57.2
## 6 chrM 1 1.57e4 87917 15542 98.7 833. 22.5 47.1
## 7 chrUn 1 1.99e7 21357254 17631961 88.7 162. 22.4 37.3
## 8 chrV 1 1.56e7 21594389 15292087 98.3 208. 22.4 57
## 9 chrVI 1 1.88e7 22998164 18593175 98.8 183. 22.5 58.5
## 10 chrV… 1 3.08e7 37101356 30240883 98.3 179. 22.4 57.9
## # ℹ 542 more rows
## # ℹ 1 more variable: samp <chr>
wg %>%
filter(!chr %in% c("chrM", "chrUn")) %>%
ggplot(aes(chr, meandepth)) +
geom_point() +
coord_flip() +
facet_wrap(~samp, nrow = 3)
# get whole genome mean depth (excluding chrUn, chrM, chrXIX, chrY)
wg_cov <- wg %>%
filter(!chr %in% c("chrUn", "chrM", "chrXIX", "chrY")) %>%
group_by(samp) %>%
summarize(weighted_mean_depth = weighted.mean(meandepth, endpos))
wg_cov
## # A tibble: 23 × 2
## samp weighted_mean_depth
## <chr> <dbl>
## 1 CH-2009 187.
## 2 CH-2010 342.
## 3 CH-2011 245.
## 4 CH-2012 355.
## 5 CH-2015 263.
## 6 CH-2016 425.
## 7 CH-2017 339.
## 8 LB-1999 107.
## 9 LB-2001 238.
## 10 LB-2002 187.
## # ℹ 13 more rows
# Does it look like even depth ratios? Would expect chrY / chrXIX depth = 0.5 3x for every Y with 50/50 so chrXIX / chrY depth of 3 -> might later years more female skewed...can I correct for this?
sex_ratio <- wg %>%
filter(chr %in% c("chrXIX", "chrY")) %>%
pivot_wider(id_cols = samp, names_from = chr, values_from = meandepth) %>%
mutate(sex_ratio= chrXIX / chrY)
sex_ratio
## # A tibble: 23 × 4
## samp chrXIX chrY sex_ratio
## <chr> <dbl> <dbl> <dbl>
## 1 CH-2009 140. 59.7 2.35
## 2 CH-2010 278. 92.5 3.00
## 3 CH-2011 196. 70.0 2.79
## 4 CH-2012 297. 86.1 3.45
## 5 CH-2015 212. 71.5 2.96
## 6 CH-2016 355. 100. 3.55
## 7 CH-2017 283. 83.0 3.41
## 8 LB-1999 92.4 24.8 3.73
## 9 LB-2001 205. 59.7 3.42
## 10 LB-2002 154. 47.8 3.22
## # ℹ 13 more rows
coverage <- left_join(wg_cov, sex_ratio)
## Joining with `by = join_by(samp)`
coverage
## # A tibble: 23 × 5
## samp weighted_mean_depth chrXIX chrY sex_ratio
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 CH-2009 187. 140. 59.7 2.35
## 2 CH-2010 342. 278. 92.5 3.00
## 3 CH-2011 245. 196. 70.0 2.79
## 4 CH-2012 355. 297. 86.1 3.45
## 5 CH-2015 263. 212. 71.5 2.96
## 6 CH-2016 425. 355. 100. 3.55
## 7 CH-2017 339. 283. 83.0 3.41
## 8 LB-1999 107. 92.4 24.8 3.73
## 9 LB-2001 238. 205. 59.7 3.42
## 10 LB-2002 187. 154. 47.8 3.22
## # ℹ 13 more rows
# Add additional depth info to rois
roi_cov <- left_join(roi, coverage) %>%
separate(samp, into = c("pop", "year"), sep="-", convert = TRUE) %>%
mutate(wg_norm_depth = meandepth / weighted_mean_depth,
chrXIX_norm_depth = meandepth / chrXIX,
Y_norm_depth = meandepth / chrY,
yrs_since_found = case_when(pop == "CH" ~ year - 2009,
pop == "SC" ~ year - 2011,
pop == "RS" ~ year - 2009,
pop == "LB" ~ year - 1985),
pop = case_when(pop == "RS" ~ "RABS",
TRUE ~ pop),
pop = fct_relevel(pop, c("RABS", "SC", "CH", "LB")),
desc = as.factor(desc),
desc = fct_relevel(desc, "NLRC5", "SYT19", "CALB2A", "HTRA1A", "MYH3C1", "MYH3C2", "MYH3C3", "MYH3C1y", "MYH3C2y", "MYH3C3y"))
## Joining with `by = join_by(samp)`
# Is it ok to use all of chrXIX for normalization?
roi_cov %>%
filter(pop == "CH" & year == "2017" & str_detect(desc, "MYH")) %>%
summarise(sum = sum(numreads),
sum_covbase = sum(covbases))
## # A tibble: 1 × 2
## sum sum_covbase
## <dbl> <dbl>
## 1 94633 65676
wg %>%
filter(samp == "CH-2017" & chr == "chrXIX")
## # A tibble: 1 × 10
## chr startpos endpos numreads covbases coverage meandepth meanbaseq meanmapq
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 chrXIX 1 2.06e7 38561406 20182037 98.1 283. 26.2 57.7
## # ℹ 1 more variable: samp <chr>
roi_cov
## # A tibble: 253 × 21
## desc region chr startpos endpos numreads covbases coverage meandepth
## <fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 NLRC5 chrXI… chrX… 2601498 2.63e6 15164 30076 95.3 67.9
## 2 MYH3C1 chrXI… chrX… 2637324 2.65e6 4770 13318 89.2 47.1
## 3 MYH3C2 chrXI… chrX… 2655758 2.67e6 3912 10519 96.8 54.3
## 4 MYH3C3 chrXI… chrX… 2667295 2.68e6 3602 10178 94.9 50.6
## 5 duplicatio… chrXI… chrX… 2665791 2.68e6 5039 16600 94.2 42.0
## 6 SYT19 chrXI… chrX… 2697955 2.73e6 11447 24674 75.5 50.9
## 7 CALB2A chrXI… chrX… 2745064 2.76e6 13620 14596 100 138.
## 8 HTRA1A chrVI… chrVI 14382418 1.44e7 25067 16913 98.7 217.
## 9 MYH3C3y chrY:… chrY 8892731 8.90e6 3282 11244 100 43.4
## 10 MYH3C2y chrY:… chrY 8904789 8.92e6 4912 10952 100 68.2
## # ℹ 243 more rows
## # ℹ 12 more variables: meanbaseq <dbl>, meanmapq <dbl>, pop <fct>, year <int>,
## # weighted_mean_depth <dbl>, chrXIX <dbl>, chrY <dbl>, sex_ratio <dbl>,
## # wg_norm_depth <dbl>, chrXIX_norm_depth <dbl>, Y_norm_depth <dbl>,
## # yrs_since_found <dbl>
If the sex ratios are 1F:1M: 3:1 chromosome ratio (chrXIX:chrY) If the sex ratios are 2F:1M: 5:1 chromosome ratio If the sex ratios are 1F:2M: 1:2 chromosome ratio
roi_cov %>%
filter(desc == "HTRA1A") %>%
ggplot(aes(yrs_since_found, sex_ratio, color = pop)) +
geom_hline(yintercept = 3, linetype = "dashed") +
geom_line(linewidth = 1.5) +
geom_point(size = 3) +
# scale_color_manual(values=c("#9970ab", "#5aae61")) +
xlab("Years since founding") +
ylab("Sex chromosome ratio (chrXIX / chrY depth)")
roi_cov %>%
filter(desc == "HTRA1A") %>%
mutate(chrXIX = chrXIX / weighted_mean_depth,
chrY = chrY / weighted_mean_depth) %>%
pivot_longer(cols = chrXIX:chrY, names_to = "chromosome", values_to="depth_1x") %>%
ggplot(aes(yrs_since_found, depth_1x, color = pop, group = interaction(chromosome, pop))) +
geom_line(aes(linetype = chromosome), linewidth = 1.5) +
geom_point(size = 3) +
# scale_color_manual(values=c("#9970ab", "#5aae61")) +
ggtitle("Sex chromsome depths") +
xlab("Years since founding") +
ylab("Normalized depth (1x genome coverage)")
roi_cov %>%
filter(chr!= "chrY") %>%
ggplot(aes(yrs_since_found, wg_norm_depth, color = pop)) +
geom_line(size = 1) +
geom_point(size = 1) +
# scale_color_manual(values=c("#9970ab", "#5aae61")) +
xlab("Years since founding") +
ylab("Normalized depth (1x genome coverage)") +
facet_wrap(~desc, nrow = 2)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
roi_cov %>%
filter(desc %in% c("NLRC5", "STY19","CALB2A", "HTRA1A", "MYH3C1", "MYH3C2", "MYH3C3")) %>%
ggplot(aes(yrs_since_found, wg_norm_depth, color = pop)) +
geom_line(linewidth = 1) +
geom_point(size = 1.5) +
# scale_color_manual(values=c("#9970ab", "#5aae61")) +
xlab("Years since founding") +
ylab("Normalized depth (1x genome coverage)") +
facet_wrap(~desc, nrow = 2)
roi_cov %>%
filter(chr!= "chrY") %>%
ggplot(aes(yrs_since_found, chrXIX_norm_depth, color = pop)) +
geom_line(size = 1) +
geom_point(size = 1.5) +
# scale_color_manual(values=c("#9970ab", "#5aae61")) +
xlab("Years since founding") +
ylab("Normalized depth (1x chrXIX coverage)") +
# ylim(0,1.7) +
facet_wrap(~desc, nrow = 2)
roi_cov %>%
filter(desc == "MYH3C3") %>%
ggplot(aes(yrs_since_found, chrXIX_norm_depth, color = pop)) +
geom_line(size = 1) +
geom_point(size = 1.5) +
# scale_color_manual(values=c("#9970ab", "#5aae61")) +
xlab("Years since founding") +
ylab("Normalized depth (1x chrXIX coverage)")
roi_cov %>%
filter(chr== "chrY") %>%
ggplot(aes(yrs_since_found, Y_norm_depth, color = pop)) +
geom_line(size = 1) +
geom_point(size = 1.5) +
# scale_color_manual(values=c("#9970ab", "#5aae61")) +
xlab("Years since founding") +
ylab("Normalized depth (1x chrY coverage)") +
# ylim(0,1.7) +
facet_wrap(~desc, nrow = 1)
## Variant analysis
af <- read_tsv(af_file) %>%
filter(pos < 2900000, pos > 2500000, `RS-2009` > 0.9) %>%
pivot_longer("RS-2009":"LB-2017", names_to = "pop", values_to = "af") %>%
separate(pop, sep = "-", into = c("pop", "year"), convert = TRUE) %>%
mutate(yrs_since_found = case_when(pop == "CH" ~ year - 2009,
pop == "SC" ~ year - 2011,
pop == "LB" ~ year - 1985,
pop == "RS" ~ year - 2009),
pop = case_when(pop == "RS" ~ "RABS",
TRUE ~ pop),
pop = fct_relevel(pop, c("RABS", "SC", "CH", "LB")),
ref_af = 1-af)
## Rows: 200098 Columns: 27
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (3): chrom, ref, alt
## dbl (24): pos, RS-2009, CH-2009, CH-2010, CH-2011, CH-2012, CH-2015, CH-2016...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
af
## # A tibble: 18,722 × 9
## chrom pos ref alt pop year af yrs_since_found ref_af
## <chr> <dbl> <chr> <chr> <fct> <int> <dbl> <dbl> <dbl>
## 1 chrXIX 2509805 C T RABS 2009 0.964 0 0.0360
## 2 chrXIX 2509805 C T CH 2009 0.978 0 0.0220
## 3 chrXIX 2509805 C T CH 2010 0.981 1 0.0190
## 4 chrXIX 2509805 C T CH 2011 0.942 2 0.0580
## 5 chrXIX 2509805 C T CH 2012 0.993 3 0.00700
## 6 chrXIX 2509805 C T CH 2015 0.88 6 0.12
## 7 chrXIX 2509805 C T CH 2016 0.907 7 0.093
## 8 chrXIX 2509805 C T CH 2017 0.832 8 0.168
## 9 chrXIX 2509805 C T SC 2011 0.994 0 0.00600
## 10 chrXIX 2509805 C T SC 2012 0.965 1 0.0350
## # ℹ 18,712 more rows
af %>%
filter(pos > 2673873, pos > 2674873, pop!="RABS") %>%
ggplot(aes(yrs_since_found, ref_af, color = pop, fill=pop)) +
stat_summary(fun.data = "mean_cl_boot", geom = "smooth", se = TRUE) +
# geom_point(data = filter(af, pop == "RABS"), size = 2.5, shape=5, alpha=0.1) +
# geom_smooth(method="lm")
# geom_line(data = filter(af, pop != "RABS"),size = 0.1, alpha=0.1) +
theme_cowplot(8) +
panel_border(color="black", size=0.75) +
theme(legend.position = c(0.8, 0.3)) +
scale_color_manual("Population", values=c("#cc79a7", "#7a669b", "#56b4e9")) +
scale_fill_manual("Population", values=c("#cc79a7", "#7a669b", "#56b4e9")) +
xlab("Years since founding") +
ylab("FW Allele frequency")
## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
af %>%
filter(pos > 2673873, pos > 2674873, pop!="RABS") %>%
ggplot(aes(yrs_since_found, ref_af, color = pop, group = interaction(pos, pop))) +
# geom_point(data = filter(af, pop == "RABS"), size = 2.5, shape=5, alpha=0.1) +
geom_line(data = filter(af, pop != "RABS"),size = 0.025, alpha=0.8) +
stat_summary(
fun = median,
geom = 'line',
aes(group=pop),
size=3,
line_type = "dashed") +
theme_cowplot(8) +
panel_border(color="black", size=0.75) +
theme(legend.position = "none") +
scale_color_manual("Population", values=c("#cc79a7", "#7a669b", "#56b4e9")) +
scale_y_continuous(limits = c(0,1), expand = expansion(mult = c(0, 0))) +
scale_x_continuous(limits = c(-.2,32.2), expand = expansion(mult = c(0, 0))) +
xlab("Years since founding") +
ylab("FW Allele frequency")
## Warning in stat_summary(fun = median, geom = "line", aes(group = pop), size =
## 3, : Ignoring unknown parameters: `line_type`
depth_plot <- roi_cov %>%
filter(desc == "MYH3C3") %>%
ggplot(aes(yrs_since_found, chrXIX_norm_depth, color = pop)) +
geom_point(data = filter(roi_cov, desc == "MYH3C3", pop == "RABS"), size = 2.5, shape=5, stroke=1.5) +
geom_line(data = filter(roi_cov, desc == "MYH3C3", pop != "RABS"),size = 1) +
geom_point(data = filter(roi_cov, desc == "MYH3C3", pop != "RABS"), size = 1.5) +
theme_cowplot(6) +
panel_border(color="black", size=0.75) +
theme(legend.position = "none") +
scale_color_manual("Population", values=c("#d73027", "#cc79a7", "#7a669b", "#56b4e9")) +
scale_y_continuous(limits = c(0,2.2), expand = expansion(mult = c(0, 0))) +
scale_x_continuous(limits = c(-.5,32.5), expand = expansion(mult = c(0, 0))) +
xlab("Years since founding") +
ylab("MYH3C3 normalized depth")
depth_plot
# Mann-Kendall test for CH read depth
mk_CH <- mk.test(pull(filter(roi_cov, pop == "CH", desc=="MYH3C3"), chrXIX_norm_depth), alternative = "greater" )
mk_CH
##
## Mann-Kendall trend test
##
## data: pull(filter(roi_cov, pop == "CH", desc == "MYH3C3"), chrXIX_norm_depth)
## z = 3.0038, n = 7, p-value = 0.001333
## alternative hypothesis: true S is greater than 0
## sample estimates:
## S varS tau
## 21.00000 44.33333 1.00000
# Mann-Kendall test for SC read depth
mk_SC <- mk.test(pull(filter(roi_cov, pop == "SC", desc=="MYH3C3"), chrXIX_norm_depth), alternative = "greater")
mk_SC
##
## Mann-Kendall trend test
##
## data: pull(filter(roi_cov, pop == "SC", desc == "MYH3C3"), chrXIX_norm_depth)
## z = 2.403, n = 7, p-value = 0.00813
## alternative hypothesis: true S is greater than 0
## sample estimates:
## S varS tau
## 17.0000000 44.3333333 0.8095238
mk_SC$p.value
## [1] 0.008130468
# Mann-Kendall test for LB read depth
mk_LB <- mk.test(pull(filter(roi_cov, pop == "LB", desc=="MYH3C3"), chrXIX_norm_depth), alternative = "greater" )
mk_LB
##
## Mann-Kendall trend test
##
## data: pull(filter(roi_cov, pop == "LB", desc == "MYH3C3"), chrXIX_norm_depth)
## z = 1.6083, n = 8, p-value = 0.05388
## alternative hypothesis: true S is greater than 0
## sample estimates:
## S varS tau
## 14.00000 65.33333 0.50000
# Variant frequency plot
af_plot <- af %>%
ggplot(aes(yrs_since_found, ref_af, color = pop)) +
geom_point(data = filter(af, pos == 2744769, pop == "RABS"), size = 2.5, shape=5, stroke=1.5) +
geom_line(data = filter(af, pos == 2744769, pop != "RABS"),size = 1) +
geom_point(data = filter(af, pos == 2744769, pop != "RABS"), size = 1.5) +
theme_cowplot(8) +
panel_border(color="black", size=0.75) +
theme(legend.position = c(0.8, 0.2)) +
scale_color_manual("Population", values=c("#d73027", "#cc79a7", "#7a669b", "#56b4e9")) +
scale_y_continuous(limits = c(0,1), expand = expansion(mult = c(0, 0))) +
scale_x_continuous(limits = c(-.5,32.5), expand = expansion(mult = c(0, 0))) +
xlab("Years since founding") +
ylab("FW allele frequency")
af_plot
# Save allele frequency plot
ggsave(
out_file_af,
plot = af_plot,
scale = 1,
width = 3,
height = 2.5,
units = c("in"),
dpi = 300,
)
final_plot1 <- plot_grid(depth_plot, af_plot, ncol = 1, rel_heights = c(1, 1), labels = c("B","C"), label_size = 12)
final_plot1
final_plot2 <- plot_grid(NULL, depth_plot, nrow = 1, rel_widths = c(1, 1), labels = c("E","F"), label_size = 12)
final_plot2
ggsave(
out_file,
plot = final_plot2,
scale = 1,
width = 6.5,
height = 1.75,
units = c("in"),
dpi = 300,
)
final_plot_h <- plot_grid(depth_plot, af_plot, nrow = 1, rel_widths = c(1, 1), labels = c("B","C"), label_size = 12)
final_plot_h